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Abstract 

Magnetic flux penetration in superconductors involves a rich variety 
of subtle phenomena, much of which is still poorly understood. Here 
these complexities are studied by formulating the Ginzburg-Landau 
equations as a lattice gauge theory. Their solutions are compared and 
contrasted with the (heuristic) Landau model of type I superconduc- 
tivity, and the (perturbative) Abrikosov model for type II supercon- 
ductors. Novelties arise as the continuum limit is approached, related 
to an effect discovered by Hofstadter. Various cautionary remarks per- 
tinent to large-scale simulations are made. 
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1. Prolegomena. 

Accurate prediction of magnetic flux penetration patterns in supercon- 
ductors poses a formidable challenge. The pictures observed are the result of 
competition between several configurations which are degenerate in energy 
(or very nearly so). Intricate structures on many length scales can result. 
Moreover, the mechanics of multivortex systems have implications for ele- 
mentary particle theory and cosmology, so the problem is of fairly general 
interest. 

The underlying complexity of the situation mandates a nonperturba- 
tive approach. Here the Ginzburg-Landau equations are recast as a lat- 
tice gauge theory, and magnetic flux penetration patterns are determined. 
Comparisons are made with the classic early models of Landau (for type I 
superconductors) and Abrikosov (type II superconductors). Discontinuities 
like those in the Hofstadter "butterfly" pattern arise while approaching the 
continuum limit and are discussed in depth. 

2. Historical background and synopsis. 

Magnetic flux penetration in superconductors has been studied for over 
half a century. It is thus appropriate to recapitulate some of the history 
of the problem as its mathematical structure is set down. Most of the 
results relevant to this paper stem from the efforts of three authors: Landau, 
Abrikosov, and Hofstadter. Their contributions are sketched as the problem 
unfolds, 
a. Landau 

The geometry of the problem considered here is the same as the one used 
originally t 1 ' 2 ' by Landau in 1937. A square flat plate of superconductor is 
situated in the xy plane, and a perpendicular field H lies in the z direction. 
The magnetic field approaches a constant at large z. 

This early work of Landau predated the Ginzburg-Landau equations 
' 3 1 by 13 years. It missed the essential difference W between type I and 
type II superconductors and ignores important effects of flux quantization. 
Still, it survives as a textbook ^ model of the intermediate state in type I 
superconductors. 

Landau predicted t 1 ' 2 ' that the magnetic field must penetrate a square 
plate of what is now called a "type I superconductor" in a pattern of stripes. 
In the above geometry, the observed pattern in the xy plane is independent 
of x and periodic in y (or vice- versa), requiring a complete spontaneous 
breakdown of the discrete [x<-> y] symmetry. Subsequent experimental work 
I 5 1 revealed that the true situation in type I superconductors is far more 
baroque, involving patterns on many length scales. The patterns do, how- 
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ever, typically possess a degree of elongation, in qualitative agreement with 
Landau. Indeed, if the magnetic field is applied at an oblique angle, the 
domain patterns align ^ in a fashion resembling his model. 

A better theoretical understanding of superconductors follows from the 
Ginzburg-Landau formalism. The major assumption of the approach ^ is 
that the free energy density f(x) of a superconductor has an expansion 

/(x) = f n + i|D^(x)| 2 + ^[(IV(x)l 2 - l) 2 - 1] + ^ 2 H 2 (x) (1) 

where tp (x) is the order parameter (i.e., the superconductor wave function), 
H(x) is the magnetic field, D = d + iA is the covariant derivative and f n is 
the free energy of the normal state. Units are chosen so as to measure 

x in units of £ = (^o/27rB c2 )^ 
A in units of ^B C 2 
/ (x) in units of 5 2 2 /27TK 2 

H(x) in units of B C 2 (2) 

where £ is the temperature-dependent coherence length and <I>o = 27rfoc/e* 
(e* = 2e) is the elementary flux quantum. Minimization of the free energy 
gives its Euler equations of motion: 

D 2 V> + V- MV = (3a) 
+ K 2 [d 2 A — d{d • A)] = J (3b) 

J = im[V>*DV>] (3c) 

Note the appearance of the Abrikosov parameter n. This Ginzburg- 
Landau formalism is quite general ^ and can be derived ^ from the BCS 
theory. 

In the present work, the Landau geometry is used exclusively. The order 
parameter ip and vector potential A depend only on x and y in the deep 
interior of a large thick superconducting plate. By use of this geometry and 
current conservation, the Euler equations simplify to 
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d 2 r - -3 + r - r 3 = (4a) 
e ab da(Jb/r 2 ) = F (4b) 

Jb = ~K 2 e bc d c F (4c) 

where e ab is the familiar Levi-Civita antisymmetric tensor (note that J a can 
be eliminated from these equations). Here, 

r(x,y) = \ip(x,y)\% (5a) 



F(x, y) = e ab d a A b (x, y) (5b) 

(both are real scalars) and all indices are two-dimensional. The microscopic 
flux density B is in the z direction, 

B(x,y) = e z F(x,y) (5c) 

b. Abrikosov 

The parameter k determines whether a superconductor is type I (k 2 < ^) 
or type II (k 2 > i). At the boundary point (k 2 = i), a partial integral of 
Eqs. (4) exists, sometimes called the "Sarma" ^ or "self-dual" 1 10 1 solution: 

- d 2 In r = 1 - r 2 = F (6) 

The existence of these two types of superconductor was first postulated 
by Abrikosov W. He solved the Ginzburg-Landau equations perturbatively, 
taking r and F to be periodic on a rectangle of size (Ax, Ay) coherence 
lengths. Each rectangle is penetrated by v units of flux 



Ax • Ay ■ B ■ f = z^$ (7) 

where v is an integer and B is the xy spatial average of the microscopic flux. 
The Abrikosov solutions are the leading terms in an expansion in (B c i — B) 
about the linear limit, with B C 2 the critical field above which the material 
goes normal (note that B C 2 = 1 here). To second order in (1 — B), the free 
energy eq. (1) is 
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f(x, y)-f n = \k 2 B 2 - \k\1 - Bf/[l + (2k 2 - + • • • (8) 
where the bar denotes xy spatial average and 

B = F 

F(x,y) = F - -n 2 [r 2 (x, y) - r 2 ] H (9) 

(compare eqs. [6]). 

For type II superconductors (k 2 > ^), the perturbative free energy eq. 
(8) is minimized when (3 is smallest, 

Type II :/3 = /3 min « 1.1596 (10) 

leading to the prediction of a triangular lattice of flux tubes W. In the 
type I case, (k 2 < ^), the perturbative formula eq. (8) predicts I 11-14 ! a 
complicated series of patterns with large (3, implying large fluctuations: 

Type I : » 1/(1 -2k 2 ) > 1 (11) 

Although these type I patterns are generally elongated, they differ in detail 
from Landau's crude model (see 2a). The point is that the superconductor 
flux density r 2 behaves like a quantum "phase space" distribution I 11-12 !. 
The x and y coordinates of the plate are essentially Fourier conjugates, like 
position and momentum in quantum mechanics. The "uncertainty princi- 
ple" underlying eq. (7) implies that a flux distribution which is independent 
of one coordinate must be sharply localized in the other, rather than the 
periodic function envisaged by Landau. Landau's model of type I super- 
conductivity is thus (oddly enough) inconsistent with the Ginzburg-Landau 
equations. 

It is obviously imperative to go beyond the perturbative formula eq. (8). 
Yet a numerical simulation must confront a novel phenomenon first reported 
by Hofstadter I 15-17 ! and elucidated in the next section, 
c. Hofstadter 

The Hofstadter phenomenon becomes relevant when the Ginzburg-Landau 
equations are formulated on a lattice. Define complex fields ip(m, n) and real 
fields F(m, n) on sites (m, n) of a lattice (x = ma, y = na; m, n = 1, ■ ■ • , L) 
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with lattice spacing a. When the partial derivatives in eqs. (3a) and (4c) 
are replaced by covariant differences, the resulting lattice equations are 

ip{m + 1, n) + tp(m — 1, n) + U(m, n)ip(m, n + 1) 
+U*(m, n — l)ip(m, n — 1) 
= [e-a 2 |V>(m,n)| 2 ]Y>(m,n) (12) 

K 2 [F(m + 1, n) — F(m, n)\ = J y (m, n) 

= Im[ip*(m, n)U(m, n)^(m, n + 1)] (13a) 

K 2 [F(m, n + 1) — F(m, n)\ = —J x (m,n) 

= —Im[ip*(m,n)ip(m + l,n)] (13b) 

U(m,n) = U(m— l,n) x exp[ia 2 F(m, n)] (13c) 

where e = 4 — a 2 and the gauge is chosen so that the vector potential lies in 
the y direction. The lattice spacing a is measured in units of the continuum 
coherence length. Since the desired solutions fill the xy plane, F and \ i[>\ 2 are 
taken periodic on a square of size LxL lattice spacings. Then the periodicity 
of the physical currents J(m, n) implies that ip is in general quasi-periodic 
rather than periodic: 

ip(m, L) = 7p(m, 1) 
ip(L, n) = ip(l, n) 

L n 

x exp[-ia 2 ^ F {m! ,n')} (14) 

m'=l n'=l 

Since ip(m + L,n + L) is single- valued, 

L L 

a 2 ^^ F(m, n) = 2tt P = a 2 L 2 B (15) 

m=l n=l 
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so B is the flux density per elementary plaquette in units of the continuum 
B c2 . Here the integer p gives the number of flux quanta penetrating the 
L x L large square. 

Equations (12-15), though simple in appearance, imply a plethora of 
strange phenomena. This complexity can be illustrated by a calculation of 
the critical magnetic field B C 2 on a finite lattice. When B exceeds this value, 
the material goes normal, and the only sensible solution to these equations 
is the trivial one where ip(m, n) is vanishes. Near this limit, ip(m, n) is small, 
and the nonlinear terms can be neglected. Then 



Fo(m,n) = constant = B = 2irp/a 2 L 2 

P 

Uo(m,n) = exp[27rij-^m] (16) 
The equation for ip(m,n) can be simplified by separating variables, i.e., 

ipo(m, n) = V gi(m) exp[ ] (17) 

1=1 L 

When the nonlinear terms in eq. (12) are dropped, it becomes 

gi{m + 1) + gi(m - 1) + 2 cos[2ir(IL - pm)/L 2 ]g I (m) = e gi(m) (18) 

which is Harper's equation [18]. The boundary condition eq. (14) implies 
that 

giim + L) = gi- p (m) (19a) 

so that 

g I (m + L 2 /p)=g I (m) (19b) 

(it may be necessary to define a "superlattice" [15] if L 2 jp is nonintegral) . 
Thus the natural periodicity of the giira) is 

L 2 /p = l/a (20) 

where a is the (rational) number of flux quanta per plaquette. 

The lattice bulk critical field B c2 (a) is determined from the largest eigen- 
value e max {a) of eq. (18) via eq. (16): 
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B c2 (a) = 27ra/[A-e max (a)} (21) 

Note that B C 2{a) is a function of a alone. 

The Hofstadter phenomenon [15] occurs when the eigenvalue spectrum 
of eq. (18) is calculated. The result for e(a) is a very striking discontinuous 
"butterfly" pattern, with an intricately organized hierarchical fine structure. 
From eq. 15, 

a 2 = (f )a (22) 

The continuum limit occurs when the lattice spacing approaches zero at 
fixed B, implying that the limit of interest is a — > 0. But the expected 
continuum limit 

lim B c2 (a) = 1 (23) 

is not obtained smoothly, instead occurring along the discontinuous upper 
boundary of the Hofstadter butterfly. Needless to say, this is unsettling 
behavior for the continuum limit of a lattice gauge theory, for presumably 
all thermodynamic functions (and not just B c2 ) will display rough structure 
as the continuum is approached. 

Two useful properties which follow from eq. (18) are 

e ma x{a) = e max (l - a) (24a) 



emax (a + N) = e max (a) (24b) 
where N is an arbitrary integer. Then 

01 

B c2 {a) = B c2 (l - a) (25) 

1 — a 

From ref. [17] values of e max (a) can be extracted (see Table 1). A plot 
of B c2 (a) versus a is given in Fig. 1. Note that away from a = 1/2. 

B c2 (a) w 1/(1 -a) 
e ma x(a) w 4-27ra(l-a) (26) 
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(though neither is ever a continuous function) and B C 2 (a) increases without 
bound as a approaches one. When a equals one, eqs. (12) have the trivial 
solution 



tp(m,n) = 1 
F(m,n) = B 



(27) 



[essentially equivalent to having no magnetic field, viz. eq. (24b), as 



3. Numerical solution of the lattice equations. 

a. Method 

The lattice Ginzburg-Landau equations, eqs. (12) and (13), are readily 
accessible to numerical simulation [19,20], though it is well to remember the 
cautionary remarks of section 2c. The scheme used here is particularly sim- 
ple. First, an initial choice of the U(m, n) and ip(m, n) is made for a given 
average B. Then eq. (12) is solved by relaxation [i.e., each of the if)(m,n) 
is determined from its nearest neighbors]. One hundred sweeps through the 
lattice prove sufficient. Given the old U(m, n) and ifj(m, n), the new F(m, n) 
are determined from eq. (13) via a finite Fourier transformation in a single 
step with B, the average of F(m,n), held fixed. The new U(m,n) are de- 
termined, and the process is repeated. Typically about 2,000 loops through 
the whole algorithm suffice. The usual checks using different starting con- 
ditions were made. An easy and adequate initial condition is F(m,n) = B; 
ip(m, n) = 1. 

The limit of interest is that of small a, by eq. (22). Yet from eq. (20), 
the natural periodicity of the system is 1/a. Thus, 1/a = L 2 /p > L. The 
optimal choice for p is therefore p = L = 1/a, and it is used here unless 
otherwise noted. 

b. Vortex arrays in type II superconductors. 

One characteristic signature of type II superconductors is the triangular 
lattice of flux tubes predicted by Abrikosov. The p maxima of F(m, n) per L 
x L periodic square are easily seen; their patterns are displayed in fig. 2 for 
parameter values B = 0.9, k 2 = 10 and various a. Since a triangular lattice 
involves irrational tangents, it can never fit exactly on a square lattice; yet 
the arrays in fig. 2 form fair approximations to a triangular lattice. In fig. 
3 the squared distances d 2 between lattice points (taken in units of squared 
lattice spacing) are plotted versus 1/a along with the Abrikosov value d 2 = 



4/3/a. Reasonable agreement is obtained, though scatter is large. The 



B/B c2 {l) = 0]. 
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lattice spacings are equal in the x and y directions, so these triangular 
"Abrikosov" lattices are a true nonperturbative prediction (compare [20]). 

For these same parameter values the xy average \ip\ 2 = p is plotted versus 
1/a in fig. 4. As the a — ► continuum limit is approached, the value of 
p scatters discontinuously [as did B c2 (a), cf. fig. 1] but clearly approaches 
a limit. Comparison with continuum values thus requires fairly small a for 
respectable results. The values for [3(a) [cf. eq. (10)] are much better: e.g. 
0(1/24) = 1.1596. 

It is well to note the existence of defect structures in lattice patterns. 
Recall from eq. (20) that the natural periodicity of a pattern is 1/a lattice 
spacings. If the system size L is not an integral multiple of 1/a, defect 
structures due to the period mismatch can occur [fig. 5(a)]. Incomplete 
equilibration can also produce defects [fig. 5(b)], in this case a superimposed 
triangular lattice of hexagonal defects. [Other parameter values are the same 
as in fig. 4]. The first problem can be eliminated and the second reduced 
by using p = L = 1/a. 
c. B versus H. 

The difference between type I and type II superconductors can be high- 
lighted by comparing B, the magnetic flux density inside the superconduc- 
tor, with H, the applied magnetic field. Here B is an input parameter given 
by the spatial average of F(m, n), while H in the present units is defined by 

[cf. eq. (8)] and can be calculated with an elegant virial theorem [21]. It is 
important to note that the assumption that F(m, n) is periodic on an L x L 
cell implies constraints by eqs. (13), 

L L 

J2Jy(™,n) = J2 J *( m > n ) = ° ( 29 ) 

m=l n=l 

For small enough B, eqs. (29) are typically violated, implying that in this 
limit the only nontrivial solutions to the Ginzburg-Landau equations involve 
widely-separated magnetic vortices. Thus the plots given in fig. 6 do not 
continue to B = 0. 

In fig. (6a), the type II case k 2 = 10 is shown. Note that, as expected, 
H is larger than B and extrapolates to a finite value H c \ as B tends to zero. 
At the boundary point k 2 = \ between type I and type II superconductors, 
[fig. (6b)]: 
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H(B) = B, B > B c2 

= H c2 ,B<B c2 (30) 

while the flux patterns generated are the same as in the type II limit. 

The B(H) curve for type I superconductors is more subtle. Figure (6c) 
displays this function for k 2 = 0.35. Note that B(H) is double- valued, with 
turning points at H = H c2 and H = H c > H c2 . As is well-known, type I 
superconductors exhibit "supercooling," leading to a hysteresis loop in the 
physically realized B{H). Bulk superconductors become normal at H = H c , 
while the normal state does not become superconducting again until H is 
lowered to H = H c2 < H c . The curve fig. 6(c) is metastable for B < B c2 , 

Metastability also occurs in the "effective potential" of quantum field theory 
and in the equation of state of thermodynamic systems (see, e.g., ref. [22]). 
There, as here, the proper behavior of the system can be determined by a 
"Landau construction," shown as dotted lines in fig. 6(c). 

The size of the metastable region is determined by the ratio H c /H c2 , 
which is 1/(s/2k) in the continuum. The metastable region thus becomes 
more pronounced as k decreases. The Ginzburg-Landau formalism for type 
I superconductors [11-14] is probably only valid for H near H c2 , and true 
experimental predictions for equilibrium structures are likely best obtained 
with "simulated annealing" methods such as that applied in [20] to type II 
superconductors. Flux patterns extracted from the lattice Ginzburg-Landau 
equations in the metastable region are generally irregular elongated nucle- 
ation lumps whose periodicity is that of the L x L system. Thus the thickness 
of the superconducting plate provides an important scale [11]. Although for 
reasons discussed above [in sec. 2(a)] Landau's model of flux penetration in 
type I superconductors is incorrect in detail, it may yet have some qualitative 
validity. 



4. Conclusions . 

As should be evident, the study of magnetic flux penetration in super- 
conductors is a fascinating and difficult subject. Here the phenomenon was 
studied by formulating the Ginzburg-Landau equations as a lattice gauge 
theory, following a review of theoretical expectations based upon heuristic 
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reasoning (Landau) and perturbation theory (Abrikosov). Taking the con- 
tinuum limit of the lattice gauge theory was a more subtle operation than 
expected, involving novelties first discussed in depth by Hofstadter. 

Nevertheless, results familiar from continuum theory were readily ob- 
tained for both type I and type II superconductors. In the latter case the 
expected triangular "Abrikosov" lattice of flux tubes was obtained without 
the traditional recourse to perturbation theory. 

The "intermediate state" of type I superconductors proved to be most 
remarkable. Here the solutions of the classical Ginzburg-Landau equations 
of motion were shown to predict a region of metastability, which perforce 
limits their domain of validity. Flux penetration in type I superconduc- 
tors seems to be controlled more by the physics of metastability than by 
the Ginzburg-Landau paradigm, though predictions can be made [11,12]. 
Time-dependent theoretical [23] and experimental [24] studies should prove 
important here, as well as studies using "simulated annealing" methods [20]. 

Possible practical applications of sensitive metastable phenomena in the 
superconducting intermediate state include a new type of "dark matter" de- 
tector [11,12] for particle physics. It may also be useful to view intermediate 
state superconductors as a giant array of Josephson junctions with dynamic 
boundaries. 

Supercomputer time provided by the Department of Energy proved use- 
ful at several stages of this work. This manuscript was typeset in LaTeX by 
Toni Weil. 
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Table 1. e max (a) for selected a. 



QL ^max (oQ 

1/2 2^/2 

1/3 1 + ^3 

1/4 2^2 

1/5 2.96645 

1/6 (5 + V21) 1/2 

1/8 [6 + (12 + 8 ^2) 1/2 ] 1/2 
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Figure Captions 



Figure 1 

Critical field B C 2(a) versus a (points). The line is the function 1/(1 — a), 
displayed for comparison. 

Figure 2 

Unit cells of vortex lattices for k 2 = 10 and a = 1/8, 1/12, 1/16, 1/20, 
and 1/24; respectively. Lines are drawn to guide the eye. 

Figure 3 

Squared distance (in units of a 2 ) d 2 between vortices plotted versus a for 
regular lattices. Dots show d 2 versus 1/a (N.B. a = 10 and 13 correspond 
to square lattices). Abrikosov result ^4/3/a shown as straight line. 

Figure 4 

Average \ip\ 2 = p versus 1/a for k 2 = 10 and B = 0.90. 
Figure 5 

Vortex patterns with defects (shaded). 

a) q = 1/8, L = 12 (periodicity mismatch) 

b) a = 1/8, L = 16 (insufficient equilibration). 
Lines between vortices are drawn to guide the eye. 

Figure 6 

B versus H for a = 1/12 

a) k 2 = 10 

b) k 2 = 0.5 

c) k 2 = 0.35. 

Solid lines denote the calculated curve; dashed lines the "Landau construc- 
tion." 
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